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In this work, we examine the effect of bulk viscosity on elliptic flow taking into account the critical 
behavior of the EoS and transport coefficients near the QCD phase transition. We found that the 
Pt dependence of V2 is quantitatively changed by the presence of the QCD phase transition. Within 
reasonable values of the transport coefficients, vi decreases by a factor of 15% at low px (< 1 
GeV). However, for larger values of pr (> 2 GeV), the interplay between the velocity of sound and 
transport coefficient near the QCD phase transition enhances V2- We further point out that Grad's 
14 moments approximation cannot be applied for the calculation of the one-particle distribution 
function at the freeze-out. 
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I. INTRODUCTION 



Among many novel discoveries in the studies of relativistic heavy-ion physics, one of the remarkable facts is the 
presence of strong collective flows, which are well described using the ideal hydrodynamic model [l|. In this sense, 
hydrodynamics is expected to be an important tool to investigate the properties of the QCD matter. If such a scenario 
is established, we may further use the hydrodynamic model to infer the initial condition just after the collision from the 
final state observables. However, there are still several open questions in the hydrodynamic description of relativistic 
heavy ion collisions [|J and more extensive studies are necessary. The effect of dissipation is one of such problems. 

The effect of shear viscosity in heavy ion collisions has been investigated by several authors [1, 0, H, IE 0> IS HI , 
and it was found that the elliptic flow parameter V2 can be significantly affected. On the other hand, the effect of 
bulk viscosity has not been analyzed in detail yet. Roughly speaking, while the shear viscosity acts as the resistance 
against the deformation of a fluid element, the bulk viscosity acts against the expansion or compression of the fluid. 
Thus, one may naturally expect that shear viscosity affects directly V2, which characterizes the spatial anisotropy of 
the dynamics. However, this does not mean that the bulk viscosity plays a secondary role. Because of the anisotropic 
expansion and compression of the produced matter, the bulk viscosity should alter the anisotropic flow as well. 

Another important factor is the behavior of the bulk viscosity coefficient, £, near the QCD phase transition temper- 
ature T c . Recently, lattice QCD simulations [hIEI suggest that £ is strongly enhanced near T c , although the shear 
viscosity coefficient is suppressed [ljj. A similar enhancement will appear even below T c , as is shown in where 
( in the hadronic phase increases also quickly towards T c , Considering the fact that the fluid dynamics slows down 
near T c due to the small sound velocity c s , the critical behavior of £ discussed above may play a crucial role on vi . 

As far as we know, the effect of the bulk viscosity on v<2 was calculated in [l5| where the effect of the phase 
transition is also considered. As a matter of fact, a quantitative estimate of the effect of the phase transition in 
dissipative hydrodynamics is not straightforward because the above critical behavior of the bulk viscosity should be 
consistently constructed with the equation of state (EoS). 

As a reliable EoS, the first choice would be the results of the lattice QCD calculation for the QGP phase and of 
the hadron resonance gas for the hadron phase. However, the connection of these two EoS's is not trivial because no 
reliable information is available near T c . In particular, c s near T c is very sensitive to the way that the two phases are 
connected and it is important to check how these ambiguities influence vi- 

The purpose of this paper is to investigate the effect of bulk viscosity on the elliptic flow within a 2+1 dimensional 
hydrodynamic modeling, incorporating the critical behaviors from c s and Q near T c . As was pointed out already 
H [H v 2 is affected also by the dissipative corrections on the one-particle distribution function at the freeze-out. 
However, we will not consider these corrections in this work because the 14 moments approximation, which is the 
method commonly employed to evaluate this correction, is inapplicable for the freeze-out values of bulk viscosity, as 
we will discuss later. 

This paper is organized as follows. In Section|TTl we briefly review the equations of causal dissipative hydrodynamics 
in the hyperbolic coordinate system and present the SPH (smoothed particle hydrodynamics) parameterization for the 
numerical implementation. In Section IIII1 we set up the parameters of our calculation, such as the initial condition, 
the EoS and the transport coefficients. We futher discuss the freeze-out process and show the intrinsic problem of 
the 14 moments approximation. The numerical results for the effects of bulk viscosity on v-i are given in Section [IVI 
Concluding remarks are given in Section |Vj 
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II. RELATIVISTIC DISSIPATIVE HYDRODYNAMICS 
A. Memory Function Method in a Hyperbolic Coordinate System 



As is well known, the application of dissipative hydrodynamics in the relativistic regime requires some caution. 
Recently, it became clear that the instability of the hydrodynamic model can be induced by the violation of relativis- 
tic causality Thus, the relativistic Navier-Stokes (Landau-Lifshitz) theory [I?) is essentially unstable and not 
applicable for the complex systems generated in the relativistic heavy ion collisions. Several formulations of relativistic 
dissipative hydrodynamics have been proposed so far to solve the problem of acausality and instability (20j |. 

In this work, we will use the memory function method fl8l . fl9| . This theory is formulated in such a way that 
the magnitude of the bulk viscosity has naturally a lower bound and we can carry out stable numerical simulations 
even for ultra-relativistic initial conditions. Here, we present the basic equations of this formalism in the hyperbolic 
coordinate system. 

For simplicity, we consider the case of vanishing baryon chemical potential where only the conservation of the energy 
and momentum is required. For a general metric g llUl we have 



1 



-A. 



r=~gT^) + I^T** = 0, 



(1) 



where T^ lX is the Christoffcl symbol and g is the determinant of g^. 

We use the Landau definition for the local rest frame and assume, as usual, that the thermodynamic relations 
are valid in this frame. In this work, we further ignore the shear viscosity. Then, the energy-momentum tensor is 
expressed as 



= (e + p + II) u»u v - (p + n) g 



(2) 



where, e, p, and II are the energy density, pressure, four velocity and bulk viscosity, respectively. 

As was discussed in [l9j], one of the important factors to obtain a stable causal dissipative hydrodynamics is to 
consider the deformation of the fluid element in the hydrodynamic flow. We introduce a reference density a defined 
by 

- = e = u% 

where 9 is the expansion rate of a fluid element and ; denotes the covariante derivative. The above relation can also 
be expressed in the form of a continuity equation for a, 
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The equation for the entropy production is written as, 
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where s is the entropy density. From this, we can define the thermodynamic force F associated with the bulk viscosity 
as V~~9 (V~9 U>1 )- I n the memory function method, an irreversible current J is induced by the thermodynamic 
force as 



J(r) 
a(r) 



TO T R {t') 



7T ex P 



T dr" 



\F(t>) J(tq) 
T , r R (r")J <j(t>) a(ro) 
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T dr' 
t Tr(t')J' 



(5) 



Thus, the equation for the bulk viscosity is obtained by setting J = U and F — \f—g 1 d tl (y'—gu^), and can be 
expressed in a differential form as follows 
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(6) 



where tr is the relaxation time. 
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In hyperbolic coordinates Eqs. {l}, © an( l © are explicitly given by 



•dr 

In this case x 11 = (t, r^, rf), where 



TR7 i. n + n = - (c + r R n) (a^ - J) . (9) 



V< 2 - z 2 , 7] = itanh f^-^ ) , r T = 0,y), (10) 



and \J—g = t. The three equations above are solved in the following calculations. 



B. Smoothed Particles Hydrodynamics 



To solve numerically the relativistic hydrodynamic equations, we use the Smoothed Particle Hydrodynamic (SPH) 
method [Hj]. See [H, Eil, H2] for the application to the heavy ion dynamics. For the sake of book keeping, here we just 
write down explicitly the SPH equations in hyperbolic coordinates leaving the detailed derivation to the references 
above. 

In the SPH method, all the hydrodynamic variables are expressed in terms of a discrete set of Lagrangian coordinates 
{r a (t), a = 1, Nsph} together with a normalized kernel function W(r,h) for the smoothing procedure. The 
parameter h represents the width of W and serves as a cut-off parameter for short wavelength modes. For hyperbolic 
coordinates, this kernel function takes the form 

W(r,h) = W(r),hr,)W(r T ,h T ), (11) 

with 



Wfahjdn = 1, (12) 
Jw{v T ,h T )dv T = 1, (13) 

where the two width parameters, h n and Ht, are introduced separately for rj and rx, respectively. 

The basic procedure of the SPH scheme is to express a reference density a* — t^o and its current j =er*v, where 

7 = 1/yl — — t 2 v 2 is the Lorentz factor, as 

Nsph 

t~1o -> o*sph (*>*) = faW(r-r a (t),h), (14) 

a=l 

NsFH rlr (+\ 

J - isPH (r, i) = ^^T^ (r " Ta W ' ft)> (15) 

a— 1 

so that the continuity equation ([3]) is automatically satisfied (22J. Eqs. ([M]) and (fl~5| can be interpreted as if there 
are particles associated with each Lagrangian coordinate r Q (i), with velocity v Q (t) = dr a /dt, carrying a quantity v a 
for the reference density cr* SPH . These particles are called SPH particles. The following calculation is independent of 
the choice of v a and, hence, we set v a = 1. 

The dynamic variables in the SPH scheme arc then, 



jr Q ,u Q ,(-^ ' ' a = 1 ' "' NsPH } 



(16) 



where they represent, respectively, position, velocity, entropy, and bulk viscosity associated with the a-th SPH particle. 
The time evolution of these variables can be obtained from Eqs. ([7]), © and © by expressing them in the SPH 
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representation. Writing the entropy density s and the bulk viscosity II as 

tjs - s* SPH (r,t) = J2^(-) W(r-r a (t),h), (17) 



n ILsph = y> Q — (-) W{r-r a (t),h), (18) 
we obtain the equation for the bulk viscosity of the a-th SPH particle 

and the equation of motion Eq. ([7]) [H, [HJ 

d /(e + p + n) a \ ^ H i^ + n^ p Q + n Q I 

^spht; ( ~Ui a J = r 2^ I/jSCT* 2 + 2 <W(|r Q - ^(t)!). (20) 

\ a ' 13=1 I fog) / 

The r.h.s. of Eq. (j2"0|) is nothing but the SPH representation of the gradients of pressure and bulk viscosity. We remark 
that the l.h.s. of Eq. (|20p contains an implicit velocity dependence through the Lorentz factor 7 in the thermodynamic 
variables. After some straightforward manipulations, we get 

Mf^ = Fj, (21) 
dr 

with 

Ml = yC5™ - -g lm u lUu (22) 
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and 



Fi = Bu. + d^p + n), (23) 



dw /n \ C 

A = + ( 24 ) 



ds \T J t r 

. 7 da* 7 1 dg lm \ . n 
S = 4-7- ~ 1 + ^^«m-^— U + — , (25) 



a* dr t 27 dr J tr' 

C = e + p + Tl. (26) 

In the above equations, we calculated the four divergence of the velocity as 

_ la da* g l °Uj du 3 UjUj dg 1 ? 
' a* dr 7 dr 2j dr ' 

and, from the continuity equation, 

da* 1 . 

—jr = — ^SPH ■ vvsph - v ■ JSPH- 
dt a 

As mentioned in the introduction, we study the transverse dynamics near the central rapidity region, assuming 
the Bjorken scaling behavior in 77, that is, u v = 0. Then the dynamics is restricted on the transverse plane and is 
simplified as 

/ A \ dip 

lycfi + -u^u ] T \ -^ = B U i r ^ d t (p + n) , (27) 



with 



(28) 

\a* dr t I t r 



These are the SPH representations of our equations which will be solved. 
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III. PARAMETER SETTINGS OF THE MODEL 



A. Initial Conditions 



One of the central issues of hydrodynamical simulations of relativistic heavy-ion collisions is to infer what are the 
real initial conditions attained in the instant of the collision. Here we just assume the initial conditions suggested by 
[H> H3 for the sake of comparison. The initial energy density is parameterized as 



:„ = K { T A ( x + -,y 



&nnTb (x 



-T B \x - -,y 



B 



a NN T A (x+ 5,1/) 




(29) 



where A = B = 197 for Au-Au collisions with the nucleon-nucleon cross section a^N = 40 mb. The free parameter K 
is fixed so that the maximum energy density Eq is 30 GeV for central collisions, as is done in [f|. The term Ta (x,y) 
is the nuclear thickness function defined as 



Ta [x, y) = dz 



Po 



1 + cxp [(r - R A ) /A] ' 



(30) 



where po = 0.17 fm 3 , Ra u — 6.37 fm and A = 0.54 fm, respectively. The impact parameter is chosen as b = 7 fm. 



B. Equation of State 



We consider three different EoS's; one for an ideal gas of massless quarks and the other two corresponding to 
different connections between the QGP and hadron gas phases, as shown in Fig. [1] The dash-dotted line corresponds 
to the ideal massless quark gas, which we refer to as EoS I. The dotted line (EoS II) and solid line (EoS III) are 
both constructed by connecting the EoS from lattice QCD calculations with three flavors for the QGP phase [ll[ and 
that of the hadron resonance gas for the hadron phase. However, as was mentioned in the introduction, it is not well 
established how the lattice EoS is connected to that of the hadron resonance gas EoS. In EoS II the two phases are 
connected smoothly for a wide interpolation domain of the temperature, 100 J$ T ;$ 290 MeV, whereas in EoS III the 
connetion is done within a narrower interpolation domain, 180 ;$ T ;$ 210 MeV (see [l2j for the connection method). 
Thus, compared with the EoS II, the EoS III is more faithful, both to the lattice QCD and the hadronic resonance 
gas values except in a narrow transition region near T ~ 200 MeV. On the other hand, since the connection is rather 
abrupt, a kink-like structure appears at the QCD phase transition. 

The temperature dependence of the pressure and c s are shown in Figs. [1] and [21 respectively. One can see that, 
although the behavior of the pressure is not drastically modified, c s is still very sensitive to the choice of connection; 
c s in EoS III has a narrow minimum reaching a much smaller value than that of EoS II. The nature of the transition 
of EoS III is almost that of the first order phase transition. If the transition were of first order, c s would vanish at the 
phase transition. Since the acceleration of the fluid is directly related to c s , this difference would affect significantly 
the hydrodynamic evolution of the system near the phase transition. Thus, we perform a comparative study of the 
behavior of V2 calculated using these equations of state. 



C. Transport Coefficients 



The behavior of the transport coefficients £ and tr has not yet been established. There are several approaches to 
calculate transport coefficients. A well-known method is to use the Green-Kubo-Nakano (GKN) formula. Exactly 
speaking, the validity of the GKN formula is, however, not clear in the relativistic regime since it was developed for 
Newtonian fluids. Recently it has been shown that even in the presence of memory effects, the viscosity coefficients 
coincide with those of the GKN formula, at least, in the leading order approximation [24[. Here, we use the result 
from lattice QCD [ic| | , ca lculated with the GKN formula, for the QGP phase. For the hadron phase, we use the result 
obtained recently by [14j |. Then, the bulk viscosity coefficient shows a maximum at T c and starts to decrease almost 
exponentially as the system departs from this critical region. The coefficient vanishes in the high temperature limit 
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FIG. 1: The temperature dependence of the EoS. The dotted and solid lines represent the EoS's with the QCD phase transition. 
The dot-dashed line denotes the ideal gas of massless quarks EoS. 
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FIG. 2: The temperature dependence of c a . The dotted and solid lines represent EoS II and EoS III, respectively. The 
dot-dashed line denotes the ideal gas of massless quarks EoS. 



in the QGP phase whereas in the low temperature limit (hadronic gas), it converges to a finite value. We fitted these 
results by analytic expressions and connected them by a parabola around T c = 200 MeV. The result is 



c 



A lX 2 + A 2 x - A 3 (0.995T C >T> 1.05T C ) 

Xiexp(-(x - 1)/cti) + X 2 exp(-(x - l)/cr 2 ) + 0.001 (T > 1.057b) 
X z exp{{x - l)/cr 3 ) + X 4 exp((x - l)/<7 4 ) + 0.03 (T < 0.995T C ), 



where x = T/T c . The fitted parameters are Ai = A 3 = 0.9, A 2 = 0.25, A 4 = 0.22, a x = 0.025, a 2 = 0.13, <j 3 = 0.0025, 
(74 = 0.022, A\ = —13.77, A 2 = 27.55 and A3 = 13.45. The Fig. [3] shows this fit and compares it with the results 
from [13, [3. 

Inspired by the result for the relaxation time of the causal shear viscosity coefficient (24|, we use the following 
parameterization for the relaxation time of the bulk viscosity, 

tr = -^b n . (31) 
e+p 

In principle, the parameter &n is a function of thcrmodynamical quantities constrained by the causality condition 
[13]. Here, for simplicity, we will always consider bu as a constant. Due to this assumption, the relaxation time also 
shows critical behavior as in the case of the bulk viscosity coefficient, as is shown in Fig. [5] 
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FIG. 3: The fit for f/s. The circles and squares denote the result obtained from Lattice QCD calculation and the hadron 
resonance gas, respectively. 
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FIG. 4: The temperature dependence of the bulk viscosity relaxation time. The solid and dotted lines represent bn = 6 and 
bn = 3, respectively. 



D. Freeze out and Kinetic Corrections 



From extensive studies of the shear viscosity, it is known that there are two different origins for the dissipative 
corrections to the elliptic flow. One is the effect of viscosity on the collective evolution of the fluid itself. The other is 
the correction to the one-particle distribution function at the freeze-out hyper surface, ft was shown that the latter 
correction is more important than the former for the shear viscosity 0, Q . 

The correction to the one-particle distribution function is usually estimated using Grad's 14 moments approximation. 
However, the applicability of this approach to the heavy-ion scenario is not so obvious. In the following, we will show 
that this procedure exhibits intrinsic problems at least for the bulk viscosity. 

In the 14 moments approximation, the one-particle distribution function for bosons is expressed as 

/ = fa + 8f^f + /o(l + fo)(E Q + A>(p"«/i) + B Q (p 2 - 4(p"u p ) a ))n, (32) 

where /o = /o (p) refers to the Bose-Einstein distribution function and is the four momentum of a particle. In 
this expression, we neglected the contributions from the shear viscosity and the heat conduction. See Appendix [Al for 
details [Hj]. 

In Figs. [5] and [HI the coefficients of the viscous corrections Bo, Dq and Eq are shown as a function of m/T. For 
the sake of comparison, we plot also Bi in Fig. [5l which is the coefficient corresponding to the viscous correction 
from the shear viscosity (see Eq. (|A29j) K One can see that the bulk correction Bq is one order of magnitude larger 
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FIG. 5: The coefficients Bo and B2 as a function oim/T. See Appendix lAl for details. 
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FIG. 6: The coefficients Do and Eq as a function oim/T. See Appendix lAl for details. 



than the shear correction B 2 - By using these coefficients with m — 140 MeV and T — 130 MeV, we calculated the 
correction to the one-particle distribution function for different values of the bulk viscosity, which are shown in Fig. 
[7J In our simulations, the typical value of the bulk viscosity at the freeze-out hyper surface reaches II 10~ 2 fm -4 , 
which means that Sf is much bigger than /o itself. Then, the corrected one-particle distribution function / can be 
even negative. The negative corrections come from the normalization conditions, Eqs. (|A20j) and (|A21j) . which imply 
that 



d 3 K 



Sf = 



d K 



ESf 



0, 
0. 



(33) 
(34) 



Obviously, the 14 moments approximation is consistent only when 5f is much smaller than the equilibrium dis- 
tribution function /q. That is, the magnitude of bulk viscosity II should be, at least, smaller than 10 -4 fm -4 at 
the freeze-out for pions. This value of bulk viscosity is very small and corresponds only to ~ 0.05% of the value 
of pressure of a hadron resonance gas at 130 MeV. Usually, the procedure of freeze-out in hydrodynamics does not 
impose any constraint on the magnitude of the bulk viscosity and it is not clear that such small values of viscosity 
will be achieved. For example, in our simulations, the value of the ratio between the bulk viscosity and the pressure 
is around 10%, which seems to be a small correction but is already too large for the 14 moments approximation. This 
situation is unchanged even for heavier particles (m/T > 1). 

One may notice that Grad's method was originally formulated for a single component fluid and, to discuss the 
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FIG. 7: The correction to the one-particle distribution function for different values of the bulk viscosity II. The solid line 
represents fo and the dotted and dash-dotted lines represent 5f for II = 10~ 3 and 10~ 2 fm -4 , respectively. 



appropriate correction to the one-particle distribution function, we have to generalize it to a multiple component 
fluid. Thus, one might argue that if we consider the multi-component fluid the limitation discussed above might be 
relaxed. In this case, each hadronic species can have its own bulk viscosity and pressure, ILu) and pa), where 

n = ]Tn w , (35) 

i 

p = ( 36 ) 

i 

Here the index (i) represents the different hadronic species. In the context of the first order theory, the ratio H^/IL 
should be proportional to the ratio of the coefficients C(i)/Ci that is, 

HW = %1 = f« „ 10 -x. (37) 

n c s 

In the above expression, s is the entropy density and su) is the entropy density contribution from the i-th specie. 
Then one should use the ratio for each species, TLu)/pu), to estimate the dissipative correction. From the discussion 
above, the magnitude of the bulk viscosity for pions, should be limited to be smaller than 10 -4 fm -4 . From 

this, we conclude that H^/p^ at T = 130 MeV should be at maximum 2%, instead of 0.05%, for the 14 moment 
method to be applicable. Even though the limit of applicability is relaxed, it is still not clear that such a constraint 
will always be satisfied in a realistic hydrodynamic scenario. Furthermore, it should be noted that the generalization 
of the 14 moments method to a multiple component fluid is not trivial, neither unique. There are several different 
solutions proposed [H, [13] • 

We conclude that the applicability of the 14 moments approximation is not clear in heavy-ion collisions. For this 
reason, we will not consider this correction to the single-particle distribution function in the following discussions. 



IV. NUMERICAL RESULTS 



We calculate the elliptic flow, V2 — (cos(20)), for pions as a function of the transverse momentum with a freeze 
out temperature of 130 MeV by the Cooper- Frye method [221 ]. Unless we say otherwise, we take &n = 6. Results for 
constant values of the ratio Q/ s will be shown for the sake of comparison. 

In Fig. [51 we show V2 calculated with the ideal gas of massless quarks EoS (EoS I). Analogous to the case of the 
shear viscosity, V2 decreases as £ increases, which is expected for smooth expanding cases, since the viscosity converts 
a part of the collective kinetic energy to heat. Even when we use the (/s with the critical behavior of the QCD phase 
transition, this surppression of V2 is essentially the same, as shown by the dotted line. As a matter of fact, we found 
a direct correlation between V2 and the entropy production. In Fig. [51 the entropy production is shown for different 
values of Q/ s, and one can notice that the larger the entropy production, the bigger the suppression of V2- 
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FIG. 8: «2 as a function of pr for EoS I with different values of £/s. The solid line corresponds to the ideal result, the dashed 
and dash-dotted lines corresponds respectively to £/s = 0.08 and 0.2 and the dotted line is the result for (/s with critical 
behavior. 
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FIG. 9: The time evolution of the entropy production for different values of f/s with EoS I (in units of the initial entropy). 
The dashed and dash-dotted lines corresponds respectively to £/s = 0.08 and 0.2 and the dotted line is the result for £/ 'a with 
critical behavior. 



Next, we consider the effect of the QCD phase transition through the EoS. In Fig. [TQ1 V2 of viscous (dashed, 
dash-dotted and dotted lines) and ideal (solid line) fluids are shown for the EoS II. The dashed and dash-dotted lines 
correspond to constant £/s cases, and the dotted line represents the £/s with the QCD phase transition. Although 
the general behavior is similar to the case of EoS I, the curvatures of «2 for the corresponding cases are changed: vi is 
reduced at low px but starts to increase at high px- This effect at high px is enhanced when £/s exhibits the critical 
behavior. 

To see this effect in detail, we show V2 up to px = 3 GeV in Fig. [IT] For px > 1.5 GeV, «2 starts to increase and 
eventually surpasses the value of the ideal case. No significant difference was observed for the EoS III compared to 
the case of EoS II, although this EoS displays a sharper transition from the QGP phase to the hadron phase. As 
discussed below, this interesting feature is the result of an interplay between the viscosity and the sound velocity. 

First, note that the enhancement of i>2 at high px is not observed for the EoS I case. Thus, this effect should be 
attributed to the presence of the QCD phase transition (or cross over). In Fig. [T^] the temperature distribution 
calculated with EoS II at r = 2.6 fm is shown. The solid and dashed lines represent the ideal and viscous results, 
respectively. One can easily see that, for the viscous case, the behavior of the fluid flow is not monotone in the region 
where c s has its minimum. This is because the flow of the internal matter, which has a higher temperature, tends 
to catch up the foregoing fluid elements, generating a non-monotonous velocity field in the radial direction. This 
behavior of the velocity field causes the bulk viscosity II to become piled up, which heats the matter in this region. If 
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FIG. 10: V2 as a function of pr with EoS II. The figure shows V2 for different values of £/s, at low transverse momentum. The 
solid line correspond to the ideal result, the dashed and dash-dotted lines corresponds respectively to £/s = 0.08 and 0.2, and 
the dotted line is the result for f/s with critical behavior. 
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FIG. 11: V2 as a function of pr with EoS II. The figure shows the behaviors for different values of f/s, at high transverse 
momentum. The solid line corresponds to the ideal result, the dashed and dash-dotted lines correspond respectively to £/s = 
0.08 and 0.2 and the dotted line is the result for f/s with critical behavior. 

this happens, the gradient of II becomes dominant in the acceleration of the fluid compared to the pressure gradient 
(note that the acceleration is given by — V (p + II) ) since the gradient of the pressure is proportional to c s . Such a 
mechanism works more efficiently in the direction where the initial acceleration is large, and in consequence, in the 
direction to increase the elliptic flow. Furthermore, this effect becomes effective only when significant collective flow 
is formed near the phase transition. Thus, the recovery of «2 by this mechanism is expected for higher momentum 
particles. We consider that this is the reason for the counter-intuitive behavior of our result. Of course one should 
keep in mind that the applicability of the hydrodynamic approach becomes dubious for the description of such a high 
Pt dynamics so that such a mechanism may not necessarily be attributed directly to experimentally observed values. 

So far, all the calculations presented were implemented fixing the relaxation time parameter bu = 6. This is because 
the bu dependence for the EoS II is negligibly small. However, the bu dependence can be observed for an EoS with a 
sharper phase transition such as EoS III. In Fig. [T31 we show the elliptic flow calculated with the EoS III for different 
relaxation times, bu = 5 and 6. One can clearly see a dependence on the parameter bu at high pt > 1.6 GeV although 
the curves of vi are independent of the values of bu for small pt < 1.6GeV. 
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FIG. 12: Comparison of the temperature profiles for the ideal and viscous fluids. The viscous fluid, shown by the solid line, 
shows an accumulation of matter near the phase transition due to the effect of viscosity. 
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FIG. 13: «2 as a function of pr for EoS III. The figure shows the results for £/s with critical behavior for different values of 
bn, in the high transverse momentum region. The solid line corresponds to the ideal result and the dotted and dashed lines 
correspond to b — 6 and 5, respectively. 



V. CONCLUDING REMARKS 

In this work, we examined the effect of bulk viscosity on v 2 taking into account the critical behavior of the EoS 
and transport coefficients near the QCD phase transition. We found that the pt dependence of v 2 is quantitatively 
changed by the presence of the QCD phase transition. Within reasonable values of the transport coefficients, v 2 
decreases by a factor of 15% at low pt (< 1 GeV). However, for larger values of px (> 2 GeV), the interplay between 
c s and C/s near the QCD phase transition generates more collective flow, with even larger values than those of the 
ideal fluid. These effects can depend on the value of the transport coefficient bn when the phase transition is more 
abrupt. Of course one should keep in mind that the applicability of the hydrodynamic approach is not guaranteed 
for the high px dynamics, but this is an interesting feature and deserves further investigation. 

All the calculations of v 2 here were done without dissipative corrections to the one-particle distribution function. 
As discussed in this paper, the naive application of the usual 14 moment calculation for the bulk viscosity leads 
to unacceptable values for the corrections to the distribution function. It should be noted that the 14 moments 
approximation is not an unique choice to solve Grad's method and there are other approaches to calculate the 
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non-equilibrium distribution function [25|. It may be possible that the problem presented here can be solved using 
alternative approaches. Further investigation on this aspect is required urgently. 

Although not explored in the present paper, one interesting aspect observed in these studies should be mentioned. 
For much larger values of £/s, say > 0.5, and smaller bu (< 3), we found some peculiar behavior of v 2 at large 
pr, exhibiting large fluctuations for small changes in the transport coefficients. We identified that such fluctuations 
are related to the emergence of instabilities in the fluid evolution and that they are amplified when a sharp phase 
transition (such as a first order phase transition) is present. We also found that these instabilities cannot be controlled 
even after introducing the additional viscosity discussed in [28| . Although the shear viscosity was not considered here, 
we believe that these instabilities are an intrinsic phenomenon that appears for large bulk viscosity dynamics, such 
as the Burgers instabilities [H, H^, |3(j. If this is true, they should reflect in the observed v 2 - For instance, the 
event- by-event fluctuations in v 2 would become very large at high py. Since this effect is closely related to the nature 
of the phase transition, it would be interesting to verify this experimentally. 

We acknowledge illuminating discussions from T. Hirano and A. Monnai and thank J. Noronha for constructive com- 
ments on the manuscript. This work has been supported by CNPq, FAPERJ, CAPES, PRONEX and the Helmholtz 
International Center for FAIR within the framework of the LOEWE program (Landesoffensive zur Entwicklung 
Wissenschaftlich- Okonomischer Exzellenz) launched by the State of Hessen. 



APPENDIX A: 14 MOMENTS APPROXIMATION 

In this appendix, we discuss the Grad's method with the 14 moments approximation (3lj |. In this approximation, 
the non-equilibrium one-particle distribution function is assumed as 

/ = (exp(-y)+a)- 1 , (Al) 
y = {e + a )+K^(e tl -p u„)+K' i K 1/ e tlv , (A2) 

where the parameters e, and e^ v represent the deviation from equilibrium up to second order in momentum and a 
is —1 for bosons, 1 for fcrmions and for a Boltzmann gas. The index is applied for the quantities in equilibrium. 
Without loss of generality e Ml/ is assumed to be traceless, 

4 = 0- (A3) 

For small momentum, the non-equilibrium one-particle distribution function can be expanded around the equilib- 
rium state 

/ « /o + /o (1 - o/o) (y -y ) + O (2) , (A4) 

where 

fo = {exp(a (x 1 t)-f3o(x,t)K>*u f t(x,t)) + ay\ (A5) 
y-yo = e + JC% + 1. (A6) 

For a rarefied gas, the conserved number current and the energy- momentum tensor are expressed as 

f d 3 K 

N>i= TTTs^^fiK'T'ri' (A7) 

C d 3 K 
J (2tt) h 

Substituting Eq. (|A4p into them, we have 

m = IS + eJg + J» v e v + Jg uX e vX , (A9) 
TV = If + eJ^ + J^\x + J Ml/Ap £A P , (A10) 

where, 

J o 1 - Q " = / 7^TT;^ Q1 ■ ■ • K ° n fo (1 " tfo) • (AH) 
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The 0-th order terms are identified as the equilibrium currents, 

IS = K, (A12) 
Ig u = Tg u . (A13) 

All the remaining terms come from the dissipative corrections. The moments (|Allj) can be expanded as 

# = Jiou", 
'o 

^6 = J io 
4 V = Jmu"u v - JaiA"", 

Jq vX = J 30 u fl u l, u x - Jaxu^A^ - J 31 u x A v ^ - J^vTA 3 * 

Jg vXp = J m u»u v u x u p - J 41 u p u v A xp - J 41 u"u x A vp - J 41 u^u p A Xv - J 41 u v u x A pp 

- J 4l u v u p A pX - J 41 u x u p A pv + J 42 A»"A Xp + J 42 A pX A vp + J 42 A pp A vX , (A14) 

where 

'"-Ww/iifi^'M^ (A15) 

Jn q = (^TTyn / p^*"* (*T fo (1 - -/o) • (A16) 

In hydrodynamics, the dissipative corrections in the currents are parameterized in terms of the bulk viscosity, 
the shear viscosity and the heat conductivity. The variables e, e M and t pv can be determined in terms of these 
hydro dynamic variables with the following set of constraints 

q v = A M „iV\ (A17) 
tt m „ = P^afi (r al3 - Tf) , (A18) 

H=~A^(T^-TD, (A19) 

Ufi (N» - Ng) = 0, (A20) 
u v (T» v - T^) = 0. (A21) 

The first three relations define the irreversible currents while the last two restrictions come from the Landau choice 
for the local equilibrium reference frame. If we were using the Eckart frame, for example, we would have to use, 
instead of (jA"20)) and (|MT|) . 



(N p - N p ) = 0, (A22) 

u»u„ (T^ - T<n = 0. (A23) 

The conditions (KT7\i , (lA"T8j) and (|Al9|) imply, 

<? = -J 2l We v - 2J 31 A^u x e uX , (A24) 

^ = 2J 42 P^ Xp e Xp , (A25) 

n = eJ 2 i + hiu x e x + J 4 iu x u p e Xp - | J 42 A xp e Xp . (A26) 

A general solution can be constructed using the following ansatz, 

e = E Q n, (A27) 

e x = D ILu x + D iqx , (A28) 

£A P = B Q (A Xp - 3u x u p ) LT + Bi (u x q p + u p q x ) + B 2 ir Xp . (A29) 
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In order to know the form of the one-particle distribution function (|A4|) . we have to determine the form of the 
coefficients E , D , Di, B , Bi and B 2 . This can be done by substituting these general solutions into the equations 
discussed above, 



(J21 




+ {J31 + J4l)B 1 = 


-1, 


(A30) 






B 2 = 


1 

2 J42 ' 


(A31) 


J21E0 -+ 


JslDa + 


(3J 4 i + 5J 42 ) B = 


1, 


(A32) 




+■ J20D0 - 


h(J 3 o + hi)3B = 


0, 


(A33) 






J31D1 + J41B1 = 


0, 


(A34) 


J20E0 


+ J30-D0 - 


\- ( J40 + J41) 3-Bq = 


0. 


(A35) 



The first three equ ation s, (|A30j) . (|A31|1 and ([A32]) . come directly fro m (|A24ll . (|A25jl a nd (|A26|) . respectively. The 
last three equations, (|A33[) . (|A34[) . and (|A35[) , are consequences of Eqs. (| A20|) and (|A21[) . The solution of this set of 
equations is 

— 3Ci J21 — 3C 2 J31 — 3J41 — 5J42 ^ ^ 

Bi = j ' h \ , (A37) 
J41J21 — 

^ = (A38) 

# = -* J ? J ? ~ J ; iJ ;° = -c 2 , ( A39) 

oi>0 ^30^10 — ^20^20 

J41J21 — ^31^31 

^ = m 2 + 4 J3lJ 3 Q- J41J20 (A41) 

3i?o J30J10 — 'ha-ha 

where m is the mass of the particle. 

In this paper, we are interested only in the bulk viscosity. Then, the form of the distribution function is 

/ = fa + fa(l - af )(Eo + D (p» U J + B ( P 2 - 4( P ^) 2 ))T1, (A42) 

We remark that the 14 moments approximation was formulated for a single component gas. The generalization of 
this method for multiple components systems has not yet been done. 
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